
%% This code estimates risk premia for each permno using a K-factor model. 
%% For each factor, the model includes two jump components corresponding to large negative and positive jumps, 
%% defined over the size brackets [-inf, -0.005] and [0.005, inf], respectively. 
%% The estimates are saved in Result_K_',num2str(K),'.mat, which are later used for creating tables and figures.


clear;

addpath('../Functions')
Filepath = '../Data/Data_wide_all/'; % Sample files are included. Raw data are restricted and cannot be shared.
% These files contain the same 15-minute return data as in Data_by_permno/, 
% but are reorganized so that each file includes all firm returns for a 21-day window.

for K=[1,3,4,5,6]



H        =   2*K;
files    =   dir([Filepath,'*.mat']);
n = 292 * 567 + 351; % Total number of observations
% There are 293 files in total in Data_wide_all/: 
% the first 292 files each contain 567 rows, and the last file contains 351 rows.
qn       =   567;  % 21 day x 27 obs per day
Delta_n  =   1/252/26; % 26 intervals per day
M        =   3488; % total number of tickers
lambdat_hat = zeros(size(files,1)-1, K+H);
Vnt_hat     = zeros(K+H, K+H, size(files,1)-1);
betat_hat   = zeros(M, K+H, size(files,1)-1);
yt_hat      = zeros(M, size(files,1)-1);
R2_t        = zeros(size(files,1)-1,1);


filesT = dir(['RESULTS_K_',num2str(K),'/*.mat']); % load jump beta estimates
TT = struct2table(filesT);
TT = sortrows(TT, 'name');
filesT = table2struct(TT);

agg = [];
todarray = [];

for i = 1:size(filesT,1)
    load(['RESULTS_K_',num2str(K),'/',filesT(i).name])
    agg = [agg;beta_J_hat];
    todarray = [todarray;todj];
end

beta_J = agg;

for i = 1:size(files,1)-1
    disp(i)
    load([Filepath, num2str(1000+i), '.mat']);

    dPt_j    =   df_all(:,3:end)';
    factors = load('../Data/factors.mat');

    if K==4
    
        df = factors.final_return_matrix(:,[1:5,8,6,7,9]);

    else

        df = factors.final_return_matrix;
    end

    dFt_j    =   df(df(:,1)>=df_all(1,1)&df(:,1)<=df_all(end,1), 3:2+K)';
    Rf_j     =   df(df(:,1)>=df_all(1,1)&df(:,1)<=df_all(end,1), end)';
    Rf_j     =   ones(size(dPt_j)).* Rf_j;
    dPt_j    =   dPt_j - Rf_j;
    
    load([Filepath, num2str(1000+i+1), '.mat']);
    dPt_jp1  =   df_all(:,3:end)';
    dFt_jp1  =   df(df(:,1)>=df_all(1,1)&df(:,1)<=df_all(end,1), 3:2+K)';
    Rf_jp1   =   df(df(:,1)>=df_all(1,1)&df(:,1)<=df_all(end,1), end)';
    Rf_jp1   =   ones(size(dPt_jp1)).*Rf_jp1;
    dPt_jp1  =   dPt_jp1 - Rf_jp1;

    [lambdat_hat(i,:), Vnt_hat(:,:,i), betat_hat(:,:,i), yt_hat(:,i), R2_t(i,:)] = AJX_HF_tod(dFt_j, dPt_j, dPt_jp1,beta_J,tod,todarray,27);
end

Lambda_hat_AJX = sum(lambdat_hat,1)/(floor(n/qn)*qn*Delta_n);
Vn_AJX         = sum(Vnt_hat,3)/(floor(n/qn)*qn*Delta_n)^2;

save(['Result_K_',num2str(K),'.mat'], 'lambdat_hat', 'Lambda_hat_AJX', 'Vnt_hat', 'Vn_AJX', 'betat_hat', 'yt_hat', 'R2_t', 'qn', 'K', 'H', 'Delta_n');
% The final estimation results are provided for replicating figures and tables in the main text.
end


